% Copyright (C) 2009,2010,2011,2012  Marco Restelli
%
% This file is part of:
%   LDGH -- Local Hybridizable Discontinuous Galerkin toolkit
%
% LDGH is free software: you can redistribute it and/or modify it
% under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% LDGH is distributed in the hope that it will be useful, but WITHOUT
% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
% or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
% License for more details.
%
% You should have received a copy of the GNU General Public License
% along with LDGH. If not, see <http://www.gnu.org/licenses/>.
%
% author: Marco Restelli                   <marco.restelli@gmail.com>


function Uint = unstr_int(x,grid,base,U,varargin)
% Uint = unstr_int(x,grid,base,U)
% Uint = unstr_int(x,grid,base,U,missing)
%
% Interpolation of the FE function U on the unstructured points x
% (passed as column vectors).
%
% First a search arlgorithm is used to find which element each point
% belongs to, and then the finite element basis is evaluated. The FE
% function U must be scalar and represented as a discontinuous field
% (i.e. as a two index array).

 np = size(x,2);

 xi = zeros(grid.d,np);
 Ue = zeros(base.pk,np);
 iei = 1;
 if(nargin>4)
   for i=1:np
     [iei,xii] = locate_point(x(:,i),grid,iei,1);
     if(iei>0)
       Ue(:,i) = U(:,iei);
       xi(:,i) = xii(2:grid.d+1);
     else
       Ue(:,i) = 0;
       xi(:,i) = 1/(grid.d+1);
       Ue(1,i) = varargin{1}/ev_pol(base.p_s{1},xi(:,i));
       iei = 1;
     end
   end
 else
   for i=1:np
     [iei,xii] = locate_point(x(:,i),grid,iei);
     Ue(:,i) = U(:,iei);
     xi(:,i) = xii(2:grid.d+1);
   end
 end

 PHI = zeros(base.pk,np);
 for i=1:base.pk
   PHI(i,:) = ev_pol(base.p_s{i},xi);
 end
 
 Uint = zeros(np,1);
 for i=1:np
   Uint(i) = PHI(:,i)'*Ue(:,i);
 end

return
